Devel_int WGCNA modules expression in NT D5, D7, and D10 samples

Render this report with ~/spinal_cord_paper/scripts/Gg_devel_scWGCNA_module_analysis_render.sh.

library(Seurat)
library(WGCNA)
library(tidyr)
library(ggplot2)
library(stringr)
library(patchwork)
library(tidyverse)
library(cowplot)
library(pheatmap)
library(gridExtra)
source("~/Neuraltube/scripts/heatmap4.R")

Load individual seurat and test WGCNA data

The individual data sets are the Day 5 (Gg_D05_ctrl_seurat_070323), Day 7 (Gg_D07_ctrl_seurat_070323), and Day 10 (Gg_ctrl_1_seurat_070323) chicken spinal cord sets. The test WGCNA data are the modules calculated on the integrated data set of all three stages.

se_path <- c("Gg_D05_ctrl_seurat_070323",
             "Gg_D07_ctrl_seurat_070323",
             "Gg_ctrl_1_seurat_070323")

clust_col <- read.csv("~/spinal_cord_paper/annotations/broad_cluster_marker_colors.csv") %>% 
  rename(broad = broad_cluster) %>% 
  select(-marker)

Order of the broad clusters for plotting purposes.

broad_order <- c("progenitors",
      "FP",
      "RP",
      "FP/RP",
      "neurons",
      "OPC",
      "MFOL",
      "pericytes",
      "microglia",
      "blood",
      "vasculature"
      )

my.ses <- list()
col_table <- list()
ord_levels <- list()

for (i in seq(se_path)) {
  # load the data sets
  my.se <- readRDS(paste0("~/spinal_cord_paper/data/", se_path[i], ".rds"))
  annot <- read.csv(list.files("~/spinal_cord_paper/annotations",
                               pattern = str_remove(se_path[i], "_seurat_\\d{6}"),
                               full.names = TRUE))
  
  if(length(table(annot$number)) != length(table(my.se$seurat_clusters))) {
     stop("Number of clusters must be identical!")
  }
  
  # rename for left join
  annot <- annot %>% 
    mutate(fine = paste(fine, number, sep = "_")) %>% 
    mutate(number = factor(number, levels = 1:nrow(annot))) %>% 
    rename(seurat_clusters = number) 
  
  # cluster order for vln plots
  ord_levels[[i]] <- annot$fine[order(match(annot$broad, broad_order))]
  
  # create index for color coding
  col_table[[i]] <- annot %>%
    left_join(clust_col, by = "broad") %>% 
    select(c("fine", "color"))
  
  # add cluster annotation to meta data
  my.se@meta.data <- my.se@meta.data %>% 
    rownames_to_column("rowname") %>% 
    left_join(annot, by = "seurat_clusters") %>% 
    mutate(fine = factor(fine, levels = annot$fine)) %>% 
    column_to_rownames("rowname")
  
  my.ses[[i]] <- my.se

}

names(my.ses) <- c("D05", "D07", "D10")
names(col_table) <- c("D05", "D07", "D10")
names(ord_levels) <- c("D05", "D07", "D10")

rm(my.se, annot)

# The reference WGCNA data. We can have several, if we want to test many at the same time
WGCNA_data = list()
WGCNA_data[[1]] = readRDS("~/spinal_cord_paper/output/Gg_devel_int_scWGCNA_250723.rds")
my.wsub =list()
my.wsub[[1]]= c(1:22)

# the name of each sample, as they appear in my.files and in the metadata of the combined object
my.samplenames = c("D05", "D07", "D10")

# This is just to add a little bit more sense to the modules, so that we don't get just a color. Corresponds to WGCNA_data
my.modulenames = list()
my.modulenames[[1]] = c(1:22)

Module gene correlation

Here, we do a correlation matrix / heatmap, to see which cell clusters group togheter. This helps us to make more detailed dotplots.
This part of the script can still be used to compare several WGCNA datasets in parallel.

# broad cluster color table
all_col <- do.call(rbind, col_table) %>% 
  rownames_to_column("sample") %>% 
  mutate(sample = substr(sample, 1, 3)) %>% 
  mutate(sample_celltype = paste(sample, fine, sep = "_")) %>% 
  select(c("color", "sample_celltype", "sample"))

Meta data

#Get a dataframe with annotations for all the samples and colors we need.
my.metam <- list()

for (i in seq(my.ses)) {
  my.metam[[i]] <- my.ses[[i]][[]]
  rownames(my.metam[[i]]) <- paste0(rownames(my.metam[[i]]), "_", i)
}
my.metam <- do.call(rbind, my.metam)

my.metam$orig.ident <- str_replace_all(my.metam$orig.ident, pattern =  "Gg_D05_ctrl", "D05")
my.metam$orig.ident <- str_replace_all(my.metam$orig.ident, pattern =  "Gg_D07_ctrl", "D07")
my.metam$orig.ident <- str_replace_all(my.metam$orig.ident, pattern =  "Gg_ctrl_1", "D10")

# my.metam$sample_celltype = paste0(substr(my.metam$orig.ident,7,9),"_",my.metam$seurat_clusters)
my.metam$sample_celltype = paste0(my.metam$orig.ident, "_", my.metam$fine)

my.metam <- my.metam %>% 
  tibble::rownames_to_column(var = "cell_ID") %>%
  dplyr::left_join(all_col, by = "sample_celltype") %>%
  tibble::column_to_rownames(var = "cell_ID")


# get sample colors
my.colsm = c("grey", "grey30", "black")
names(my.colsm) <- c("D05", "D07", "D10")
 

Average expression data

cell_table <- my.metam %>% 
  rownames_to_column("cell_ID") %>% 
  select("sample_celltype", "cell_ID")

avg.mod.eigengenes <- WGCNA_data[[1]]$sc.MEList$averageExpr %>% 
  rownames_to_column("cell_ID") %>% 
  left_join(cell_table, by = "cell_ID") %>% 
  column_to_rownames("cell_ID")
  

# add metadata

avg.mod.eigengenes.mean <- avg.mod.eigengenes %>%
  group_by(sample_celltype) %>%
  summarise_all("mean") %>%
  column_to_rownames("sample_celltype")

spearman correlation heatmap

annotations

# names and colors for the heatmap annotation
annot_name <- data.frame(
  "Celltypes" = all_col$sample_celltype,
  "Sample"    = all_col$sample,
  row.names = all_col$sample_celltype)
  

annot_module <- data.frame(
  "Module" = colnames(avg.mod.eigengenes)[1:22],
  row.names = colnames(avg.mod.eigengenes)[1:22]
)

pheat_col_table <- do.call(rbind, col_table) %>% 
  rownames_to_column("sample") %>% 
  mutate(sample = substr(sample, 1,3)) %>% 
  mutate(fine = paste(sample, fine, sep = "_"))

# match color table with annotation
pheat_col_table <- pheat_col_table[match(annot_name$Celltypes, pheat_col_table$fine),]

annot_col <- list(
  Celltypes = pheat_col_table$color,
  Sample = c(D05 = "#A4A4A4",
             D07 = "#515151",
             D10 = "#000000"),
  Module = str_remove(colnames(avg.mod.eigengenes)[1:22], "AE")
  )

names(annot_col[[1]]) <- annot_name$Celltypes
names(annot_col[[3]]) <- colnames(avg.mod.eigengenes)[1:22]

heatmap of module pseudobulk average expression

pdf("~/spinal_cord_paper/figures/Fig_1_devel_module_v_clusters_heatmap.pdf", width = 8, height = 10)
grid.arrange(htmp$gtable)
dev.off()

null device 
          1 

Load the integrated data set

The integrated data set on which the WGCNA is calculated. We plot it, split by the original cell types from the three original samples.

# This is a file of all the combined mouse datasets, normalized and such.
my.sec = readRDS("~/spinal_cord_paper/data/Gg_devel_int_seurat_250723.rds")

identical(rownames(my.metam), colnames(my.sec))

[1] TRUE

my.metam$sample_celltype <- factor(my.metam$sample_celltype, levels = all_col$sample_celltype)
#Set the identities of the integrated data, to the annotated clusters
my.sec = SetIdent(my.sec, value = my.metam$sample_celltype)

p1 <- DimPlot(
  my.sec,
  reduction = "tsne",
  label = TRUE,
  repel = TRUE,
  cols = all_col$color,
  split.by = "orig.ident"
  ) + 
  NoLegend()

p1

pdf("~/spinal_cord_paper/figures/Devel_split_tsne.pdf", height = 7, width = 15)

#Plot split tsne
p1

Avg. module exp. by stage

tSNE DimPlots showing the average expression of each module by stage.

AE over time

We plot the average expression of each module in the three stages and the 5 broad cell type clusters present in all 3 stages.

VlnPlots of avg. module exp. by stage and seurat cluster

colored by module

# reorder seurat clusters
for (i in seq(my.ses)) {
  my.ses[[i]]$seurat_clusters <- factor(
  my.ses[[i]]$seurat_clusters,
  levels = levels(my.ses[[i]]$seurat_clusters)[as.integer(str_extract(ord_levels[[i]], "\\d{1,2}$"))]
  )

}

vplots <- list()

for (i in seq(my.ses)) {
  
  mods <- colnames(my.ses[[i]][[]])[grep("^AE",colnames(my.ses[[i]][[]]))]
  
  vplots[[i]] <- VlnPlot(
        my.ses[[i]],
        features = mods[module_order],
        group.by = "seurat_clusters",
        stack = TRUE, flip = TRUE,
        cols = substring(mods, 3)[module_order]) +
    theme(legend.position = "none") +
    geom_hline(yintercept = 0, lty = "dashed")
  
}

vplots[[1]]
vplots[[2]]
vplots[[3]]

pdf("~/spinal_cord_paper/figures/Fig_2_AE_by_cluster_modcol.pdf", height = 20, width = 10)
vplots[[1]]
vplots[[2]]
vplots[[3]]

colored by cell type

clust_col <- read.csv("~/spinal_cord_paper/annotations/broad_cluster_marker_colors.csv")

vplots_id <- list()

for (i in seq(my.ses)) {
    
  mods <- colnames(my.ses[[i]][[]])[grep("^AE",colnames(my.ses[[i]][[]]))]
  
  vplots_id[[i]] <- VlnPlot(
    my.ses[[i]],
    features = mods[module_order],
    group.by = "seurat_clusters",
    stack = TRUE, flip = TRUE,
    fill.by = "ident",
    cols = col_table[[i]]$color[as.integer(str_extract(ord_levels[[i]], "\\d{1,2}$"))]) +
    theme(legend.position = "none") +
    geom_hline(yintercept = 0, lty = "dashed")
}

vplots_id[[1]]
vplots_id[[2]]
vplots_id[[3]]

pdf("~/spinal_cord_paper/figures/Fig_2_AE_by_cluster_clucol.pdf", height = 20, width = 10)
vplots_id[[1]]
vplots_id[[2]]
vplots_id[[3]]

MN modules

For the figures we specifically select the two MN modules and plot them as Vln plots.

VlnPlot(
    my.sec,
    features = mods[module_order],
    group.by = "sample_celltype",
    stack = TRUE, flip = TRUE,
    cols = substring(mods, 3)[module_order]) +
    theme(legend.position = "none") +
    geom_hline(yintercept = 0, lty = "dashed")



pdf("~/spinal_cord_paper/figures/Fig_2_AE_by_cluster_integrated_data.pdf", height = 20, width = 30)
VlnPlot(
    my.sec,
    features = mods[module_order],
    group.by = "sample_celltype",
    stack = TRUE, flip = TRUE,
    cols = substring(mods, 3)[module_order]) +
    theme(legend.position = "none") +
    geom_hline(yintercept = 0, lty = "dashed")

```r
# Date and time of Rendering
Sys.time()

sessionInfo()

```